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PROGRAM DESCRIPTION 


The research being performed under NASA Grant NAG 1-1361 involves a 
more clear understanding and definition of the constraints involved in the pole- 
zero placement or assignment process for multiple input, multiple output 
systems. Complete state feedback to more than a single controller under 
conditions of complete controllability and observability is redundant if pole 
placement alone is the design objective. The additional feedback gains, above 
and beyond those required for pole placement can be used for eigenvalue 
assignment or zero placement of individual closed loop transfer functions. 
Because both poles and zeros of individual closed loop transfer functions 
strongly affect the dynamic response to a pilot command input, the pole-zero 
placement problem is important. 

When fewer controllers than degrees of freedom of motion are available, 
complete design freedom is not possible, the transmission zeros constrain the 
regions of possible pole- zero placement. The effect of transmission zero 
constraints on the design possibilities, selection of transmission zeros and the 
avoidance of producing non-minimum phase transfer functions is the subject of 
the research being performed under this grant. 

PROGRESS FOR 1 JULY - 1 SEPTEMBER 1992 

EGR Activity 

In the last progress report, it was proven that the determination of 
transmission zeros was a straightforward task using an expanded definition of 
Cramer's Rule or using Gauss' algorithm for determinants. It was also shown 
that the control law that decouples the system places poles at the transmission 
zero locations. In fact, the feedback recreates the eigenstructure of the 
decoupled subsystem associated with the transmission zeros, given by equ. 1 
below: 


x-i(t) = [An • BiBa-’Aai] x i(t) + A12 x 2(t) + Bi u c (t) (1) 

Because the transmission zeros lls - An + Bi 82*^21 1 = 0 are invariant, 
the only way to reduce the exitation to x^t), if it is desired to suppress the x-i(t) 



response, is to try to decouple x-i(t) from x 2 (t) and/or to interconnect the 
controllers in such a way that the effective control effectiveness Bi is reduced, 
but this may not be desirable. 

Transmission zeros are manifest only when fewer independent 
controllers than degrees of freedom of motion exist, and the primary purpose for 
decoupling is to provide for complete design freedom of the outputs x 2 (t) of the 
system, those motions of the airplane that we wish to control with exacting 
precision. 

By decoupling, the design freedom of the output is complete and 
exacting, with many options open to the flight control system designer. The 
responses x 2 (t) = y(t) can be decoupled from each other or made to have any 
behavior desired. For instance, eigenvector assignment becomes a relatively 
straight-forward matrix algebra problem, and is shown to be a simple dynamic 
inversion approach to design. The objection to dynamic inversion (that the 
closed loop system is not robust) can also be easily overcome, at least in the 
system outputs y(t) = x 2 (t) and not just one, but an entire family of solutions is 
possible. 

Consider the part of the system that was decoupled in such a way that 
the control effectiveness matrix B 2 is nonsingular, i.e., IB 2 I * 0. 

x 2 (t) = A 22 x 2 (t) + B 2 u(t) (2) 

Let us assume that we wish to design a system that behaves exactly as 

y(t) = Ly(t) + B L u c (t) (3) 

where the matrix L can be anything, including a Jordan form L = A or a modified 
Jordan form L m , i.e. 
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or even a nonlinear matrix. In the design process to follow, the choices are 
entirely up to the designer. Because the following design technique has so 
many options, the vehicle can be made to have flying qualities taken exactly 
and precisely from the flying qualities specification MIL-F-8785(c) or MIL STD 
1797, thereby "assuring" level 1 flying qualities. The only problem is that not all 
dynamics yielding level 1 flying qualities have been defined. The design 
objective or criteria is incomplete, but that is another story for near future 
research activity. Assuming the matrix L can be specified and the desired 
control matrix B|_ is definable (more easily than L), then the most direct solution 
to this problem is to define a regulator in the error between the actual response 
x 2 (t) and the desired response y(t); i.e., 



<S(t) - (Aas + P) e(t) = 0 

(5) 

where 

e(t) = x 2 (t) • y(t) 

(6) 


e(t) = x 2 (t) - y(t) 


so the regulator is 




[* 2 (t) • y(t)] • (A 22 + P) [x a (t) - y(t)] = 0 

(7) 


and the matrix P is anything the designer wishes to specify. P need not be 
linear nor even the same dimension as A 22 . One of the configurations of P 
shown below in fact exceeds the dimension of A 22 because the integral of the 
error is included in the design. 

Substituting for x 2 (t) from equation (2) into equation (7) yields 


A 22 X 2 (t) + B 2 u(t) - y(t) - (A 22 + P)(x 2 (t) - y(t)) = 0 


( 8 ) 



and solving for the control motion u(t) will produce the required control law: 


u(t) = B 2' 1 [y(t) + Px 2 (t) - (Ass - P) y(t)] (9) 

An entire family of solutions can be obtained from equ. (9). For instance, 
by substituting for y(t) and choosing A 22 + P = U equ. (9) becomes 

u(t) = Bs-MLyO) + B L u c (t) + (L - A 22 ) x 2 (t) - A^ y(t) - (L - A 22 ) y(t)] 

= B 2 - 1 [(L - A 22 ) x 2 (t) + B l Uc(t)] (10) 

The substitution of the feedback control law of equ. (10) into equ. (2) 

yields 


X 2 (t) = A22X 2 (t) + B2B2- 1 [(L - A^) x 2 (t) + B L u c (t)] 

= L x 2 (t) + B L u c (t) (11) 


and the system behaves “exactly* as desired. If L was chosen non-interacting, 
the system would be decoupled. In block diagram form, the system is shown 
below. This architecture is often called "implicit model following." 


Uc(t) 



Fig. 1 Feedback Solution 

The decoupling feedback B^A^ has been added to the diagram to 
show the complete solution. If the plant contained as many independent 
controllers as degrees of freedom of motion (such as elevator/direct lift 
flap/throttle or elevator/thrust vectoring/throttle or even canard/thrust 
vectoring/throttle), the feedback B^A^ xi(t) would not be necessary, a complete 
dynamic inversion is possible. In general, if the matrix L of the system shown in 
fig. 1 has eigenvalues larger (i.e., farther in the Ihp) than A 22) the robustness is 
improved from that of the open loop. The feedback is normally to augment the 



aircraft natural dynamics in the "degenerative" rather than "regenerative" sense; 
stability is enhanced. 

The same criteria, i.e., y(t) = Ly(t) +B L u c (t) can also be obtained without 
using feedback at all. All that is necessary is to calculate the forces and 
moments that must be applied to the aircraft to force it to behave as defined by 
equ. (11). This is done simply by defining P = 0 in equ. (9) above. Setting P = 0 
in equation (9) yields the very simple solution 

u(t) = Ba'W) * A 22 y(t)l (12) 

indicating that y(t) and y(t) are generated independently in a computer and the 
plant controllers are driven properly by the computer outputs y(t) and y(t). The 
block diagram is shown in fig. 2 below. This architecture is often called "explicit 
model following." 



Fig. 2 Feedforward Solution 

As indicated in equ. (12) and shown in fig. 2, the configuration is entirely 
an open loop architecture, yet yields "exactly" the same result as shown in 
fig. 1, a feedback solution. The control law "gains," B 2 ' 1 and B 2 _1 A 22 are a 
function only of the plant, and these gains constitute a complete dynamic 
inversion of the decoupled part of the plant. Because this is so, the model in the 
computer can be changed at will without changing the control law in any way. 
If, for instance, the flying qualities requirements were different for each flying 
task or flight regime, changing only the model dynamics will properly change 
the dynamic response. In fact, a ground-based simulator operates in exactly the 
way described above. The only difference is that the airplane is moving, not 
standing still, so the velocity and dynamic pressure dependent derivatives A 22 
are finite and must be taken into account. 










The closed loop or feedback solution of fig. 1 and the open loop or 
"feedforward" solution of fig. 2 represent two extremes of the family of solutions 
that can be obtained. Neither system is particularly robust in the sense that 
flight condition variations of the stability and control derivatives B 2 and A 22 will 
cause deviations in the desired, response unless a lot of gain scheduling is 
used. Even with gain scheduling, the lack of low frequency robustness of the 
system shown in Fig. 2, would preclude its use in an actual airplane. However, 
robustness can be easily provided and in fact, the system shown in fig. 2, with 
feedback added, can result in a robust architecture, accounting for the variation 
of the stability and control derivatives of B 2 and A 22 . In fact, there is reason to 
believe, as shown below, that minimal gain scheduling is possible for many 
aircraft. 

The family of solutions that will yield a robust system configuration is 
given by the general solution of equ. (9). The block diagram of the system 
represented by equ. (9) (including the decoupling feedback) is shown in fig. 3 
below: 



The general solution of equ. (9), (i.e., fig. 3), shows that the matrix P acts 
on the error between the desired output y(t) and the actual output x 2 (t). If B 2 and 
A 22 are known "exactly", the error is zero and P has no function. In fact, P 
defines not only the regulation of the error between y(t) and x 2 (t) but also 
defines the perturbation response of x 2 (t) to external disturbances (gust 
alleviation). Most importantly, because P can be an integration process, low 
frequency robustness can be provided, so the response of the computer output 
and the aircraft output will not diverge. Because all the flying qualities 
requirements can be resident in the computer model, the gust alleviation or 










structural mode control function can be done using feedback without affecting 
the maneuvering flying qualities requirements. In fact, equ. 9 is simply a 
formalization of the present common practice of providing feed forward gains in 
a flight control design. 

If the computer is programmed such that the kinematics of the plant are 
resident in the computer but the aerodynamics in the computer are chosen for 
flying qualities purposes, the computer can represent a global "model following" 
system. If such a system were used for a wide flight range vehicle such as 
HSCT, the HSCT vehicle would always be at the same flight condition as the 
computer generated "model," but the HSCT vehicle would respond dynamically 
as the "level 1" flying qualities model. There is reason to believe that the 
configuration defined by fig. 3 would work well for a wide flight range vehicle. 

Because the function P acts on the error between the computer 
generated model response and the actual vehicle response, P can represent a 
robust compensator as defined by such methods as an system designed- to 
minimize the error as the stability and control derivatives A 22 and B 2 vary with 
flight condition. 

In general, any compensation network can be expanded in partial 
fraction expansion form to represent proportional, integral and derivative 
components (or designed as a PID system). To show how such a network could 
be designed for an architecture of this type, consider a PID system, resulting in 
an error control law. 


u e (t) = -Kie(t) - K 2 Je(t)dt - K3 e(t) 


(13) 



shown in fig. 4 below: 



This block diagram is drawn to highlight the effects of the PID feedback 
and the effect of this feedback on the response of the system as stability and 
control derivatives B 2 and A22 vary in flight. As shown in the figure, the accuracy 
of the design depends upon the accuracy of B 2 and A 22, the matrices involved in 
the dynamic inversion. If the feedback gains are made sufficiently large, the 
system can be made insensitive to variations in the stability and control 
derivatives, i.e., 

if K 3 »B 2 1 -> insensitive to variations in control effectiveness 

and in control surface nonlinearities 

if K-i »-B 2 ' 1 A 22 -» insensitive to variations in B ^ 1 A&, i.e., the 

stability derivatives 

Making the fair assumption that the feedback is to maintain or improve 
the stability of the vehicle, the integral of the error will guarantee x 2 (t) = y(t) in 
the long term, regulating always about the trajectory of the airplane on the 
model computer (for instance, a hypersonic NASP with ideal level 1 flying 
qualities) without the need to directly measure some air data (such as velocity) 
on the vehicle itself. Because the integral guarantees long term zero error, the 
pilot command input-response output is independent of flight condition and can 
be made constant or vary with flight condition exactly as specified (Fs vs. n z ) in 
the F.Q. specifications. Because in the long term the aircraft steady state 











response will be totally predictable, a display of the pilot command input will tell 
the pilot what the steady state flight variable changes will be regardless of the 
sluggishness of the vehicle response, resulting in a "predictor 11 display. By 
limiting commands at the model, functions such as angle of attack and sideslip 
can be limited, again useful for a NASP system. The problem solving potential 
of this type of system architecture seems significant. 

Even more advantages for this architecture can be realized. For 
instance, because the gains B2' 1 A22 represent a "division" process of 
dimensional stability and control derivatives, many of which are 
dimensionalized in exactly the same way, many of the terms of B2' 1 A22 are 

simply (constant) ratios of non-dimensional stability derivatives, such as 
Cm a /cm§ e - which can be accurately obtained in a wind tunnel. These constants 

occur along the diagonal of the matrix B2 _1 A22, so the robustness oriented gains 
would be designed to minimize the effects of off-diagonal terms, (if their effect is 
significant). Because the aircraft can be made to respond to disturbances 
entirely independently of the computer model (which doesn't respond at ail 
unless turbulence is measured and injected into the computer), the feedback 
can be tailored for gust alleviation or structural mode control or even to 
minimize wring root bending moments (to disturbances) without affecting flying 
qualities (resident in the computer). The minimization of maneuver loads can 
be dealt with in the model computer and reflected in the vehicle itself. For 
instance, if the "model" had a direct lift flap, distribution of “model" wing loads 
would be reflected in the vehicle itself. 

To the designer, perhaps the biggest advantage is the fact that most of 
the ground based simulator setup can be transferred to the actual vehicle, 
accounting only for the fact that the actual airplane is moving rather than 
standing still. Accurate evaluation of the differences in pilot opinion between 
the ground based simulator and the actual vehicle can be made. In a simulator, 
the velocity is displayed to the pilot; in the actual airplane, not only can the 
velocity from the computer model be displayed, but the airplane will be moving 
that fast. 

It is clear that the sensitive part of the type of system described herein is 
Xi (t), the dynamics decoupled from x 2 (t). This will be particularly true if the 



transmission zeros are relatively low frequency. Because the "outputs" and 
measurement set or sensors can be entirely different, it may be possible, using 
observer synthesis ideas of Bacon (Langley) to improve the robustness of the 
u(t) = -B 2' 1 A 2 i x-i(t) feedback portion of the system. These possibilities will be 
considered by EGR during the next reporting period. 

Flying Qualities Implications 

It is not to suggest that decoupling will solve all flying qualities problems. 
The consequences of decoupling can be considerable. Consider the relatively 
simple modal decoupling design approach used in the Shuttle, which used 
pitch rate feedback and proportional plus integral compensation in the error 
loop. This pitch rate pole-zero canceling scheme essentially eliminated a 
phugoid mode residue in pitch rate, but increased the phugoid contribution to 
the a(t) response. The closed loop a(t) transfer function contains a pole at the 

origin, and a step stick command will yield a smooth, well-behaved pitch rate 
response, but the angle of attack will be divergent. Since y=0-CX, the attitude 

and flight path will not be harmony, i.e., the vehicle flight path will wander if the 
pilot holds a steady attitude. A pilot wants flight path proportional to attitude, so 
this kind of system may produce a PIO problem for the pilot, particularly during 
flare and landing. 

As the astronaut pilot pulls (or if he has learned correctly, pulsed) the 
stick to flare to a new pitch angle, the lift (-Z a OC) increases because CC(t) is 

divergently increasing, and the vehicle will have a tendency to “float" down the 
runway. The increase in CC(t) increases induced drag, decreasing the velocity 
and decreasing Z a . So eventually either the vehicle will settle to the runway 
(decreasing Z a greater affect than increasing a) or the vehicle will stall. The 
pilot may have to pulse forward on the stick to get down. Pilots do not like to 
push on the stick to settle on the runway, they usually depend on the phugoid 
response to settle the airplane down, (but nonexistent in a rate command, 

attitude hold system). Also, since — = — , the pilot must cope with the 

Aq s + l/T 02 

delay or time constant (2-2.5 sec in the Shuttle) of the flight path bending 
response after initiating a pitch rate response; the astronaut pilot must develop 
a considerable precognitive predictive capability. When coupled with the 



elevon location on the Shuttle, which suggests that the pilot sits aft of the center 
of rotation (percussion), thereby producing an initial "wrong direction cue" or 
time delay, then the Shuttle PIO evident in an early Shuttle flight is no surprise; 
accuracy or "tight" flight path control is out of the question. The extensive (and 
expensive) astronaut-pilot training is justified. It would appear that a little static 
stability (i.e., a(t) feedback) would help the Shuttle settle to the runway and 
provide for a better capability to maneuver near the ground without a PIO 
tendency. 

The result is that decoupling can significantly alter the eigenvector 
configuration of the vehicle, and change response residues in such a way that 
the vehicle no longer behaves in an "airplane-like" fashion. If the kind of 
change in vehicle behavior can be produced in a relatively simple decoupling 
situation as evidenced on the Shuttle, then multivariable decoupling may be 
expected to complicate the flying qualities scenario even more. Decoupling 
may be a blessing or a curse, and much research must yet be done to define 
the blessings and the curses. 

Generally speaking, it is an objective of "reconfiguration" technology to 
switch to an alternate set of control effectors upon the failure of a control 
surface, without changing the dynamic response of the vehicle. If the original 
decoupling was obtained using primary control effectors (i.e., elevator 
pitching moment, direct left flap -> lift or flight path, and/or throttle -> velocity) the 
likelihood that the transmission zeros will be rhp is small. If a switch is made to 
a secondary set of controllers during reconfiguration to maintain the same non- 
interaction, the likelihood of rhp transmission zeros is considerably increased. It 
might be prudent to abandon the decoupling criteria during reconfiguration and 
revert to a flying qualities criteria that would result in a “conventionally flying 
airplane." 

Examples of the changes in transmission zeros as a function of controller 
sets to achieve the same decoupled behavior will be shown in the next 
reporting period. 
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The theoretical developments and underpinning of the transmission zero work are described 
in the section titled “EGR activity.” The U.B. activity includes and/or is based on this theory, 
but for brevity the reporting of these developments is not repeated here. 

The major activity of the reporting period was the creation, debugging, and transfer (to 
EGR) of detailed computer simulations to be used for testing and demonstrating the theoretical 
developments. The computer simulations represent the experimental Total In Flight Simulator 
(TIFS) aircraft, operated by the CALSPAN Flight Research Department. 

For simulation purposes, the TIFS is assumed to be flying as a rigid body in a vertical 
(longitudinal) plane, so that there are three degrees-of-freedom (horizontal translation x, vertical 
translation z, and pitch angle 6 ). Aerodynamic coefficients are consistent with Calspan reports. 
The details of the simulation follow. 


Coordinate Systems 

In developing the aircraft equations of motion, two different coordinate systems are used, 
shown in Figure 1. One is an inertial coordinate system, which is fixed to the earth and is 
considered to be a non-rotating system. This is a valid assumption since the rotation of the 
earth is negligible in most aircraft dynamic problems. The other system is fixed to the aircraft 
center of gravity and rotates along with the aircraft. This system is referred to as the body axis 
coordinate system. 


l 



Figure 1 Coordinate Axis System 


Equations of Motion 

Newton’s second law is used to derive the rigid body equations of motion, i.e., the 
conservation of both linear and angular momentum (Nelson [2]): 

= (1) 

M = -7 -H = -7- t x mv (2) 

^ dt dt^ 

The airplane is considered as a continuum of mass panicles (5m) and each elemental mass 
has a velocity (v) relative to an inertial frame. Newton’s second law is: 


1 


( 3 ) 




Figure 2 Body Axis System 


From Figure 2, the velocity of the differential mass can be expressed as: 

(4) 

Since the mass of the aircraft is assumed constant and the total mass of the aircraft is found 
by summing the mass elements, Newton’s second law may be written as: 

£>#=# = y>>» (5) 

The right hand side of Equation (5) equals zero since r is taken from the center of mass. 
Equation (5) now reduces to: 


F = 



(6) 


The moment equation is developed in terms of the relative velocity of a mass element to 
the center of mass: 

dr _ _ 

v = v c + — = v c + dj x f (7) 

at 

where u; is the angular velocity of the aircraft. Substituting this expression into the moment 
Equation (2), the total moment of momentum is: 


3 




( 8 ) 


H = ^ ( r)8m x v c 4- ^ [r x (u; x f)\8m 

Again, the left hand side of Equation (8) equals zero and the moment of momentum equation 
reduces to: 


H = [f x (w x r)]8m (9) 

A difficulty occurs if the reference is rotating, since the moment will vary with time. To 
eliminate this difficulty, a transformation of the reference frame is made to the body axis system 
(onto the aircraft). A vector A is transformed from a fixed system to a rotating coordinate 
system by: 

fixed = ^-)rot.+ lux A] ( 10 ) 

dtr j dt j l J 

The force and moment equations, transformed into the body axis system, now become: 

F = TTi-j- + m(d/ x tf c ) (11) 

at 

*~7T (12) 


Vector Components and Scalar Equations 

To obtain the time history results of Equations (11) and (12), it is necessary to express the 
vector equations into component form. The component forms of the forces, moments, gravity, 
linear and angular velocities are shown in Figure 3. 1 The corresponding equations are: 

F = F x i + F y j + F z k 
M = Li + Mj + Nk 

9 = 9xi + g y j + 9zk (13 ^ 

w = Pi + Qj + Rk 
Vp = Ui + Vj + Wk 

r = xi + yj + zk 

1 Positive sense is in the direction of the arrows 
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Equation (11) can now be expanded into scalar form: 


F x + mg x = m(u -VR + WQ^j 

Fy+mg y = m(v + UR-WP ) (14) 

F z + mg z =m(w -UQ + VP) 




Z 


Aerodynamic and Thrust Moments 



z 


Linear and Rotational (Angular) 
Velocities 


Figure 3 Vector Components (from Roskam [1]) 


After expanding Equation (9), the scalar components of the moment of momentum are: 
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( 15 ) 


H x = P^2(y 2 + z 2 )S m — Q^ j {xy)6m- ( xz)5m 

H y = — P "Y^ ( xy)6m + Q (x 2 + z 2 )6m — R ( yz)8m 
»> = -pY, (x z)6m — Q ^ ( yz)Sm 4- R (x 2 + y 2s )Sm 

These equations can be expressed in terms of mass moments of inertia about the x, y, and 
z axes: 


H x = I XX P — I xy Q — I XZ R 

Hy = —I X yP + IyyQ ~ I yZ R (16) 

Hz = -hzP — lyzQ + IzzR 

Since, for most airplanes, the X-Z plane is a plane of symmetry, the moments of inertia for 
Ixy = ly* = zero. With this assumption, and applying Equations (16) to the moment expression 
in Equation (12), the scalar moment equations can be written as: 

L — I XX P — IxzH — IxzPQ + (Izz — Iyy)HQ 

M = IyyQ 4- (/„ - I ZZ )PR + I xz {p 2 - R 2 ) (17) 

N = IzzR ~ IxzP + {lyy — Ixx)PQ + IxzQR 

where L, M and N are the moments about the X, Y and Z axes, respectively. 

Earth Fixed System and the Kinematic Equations 

The force and moment equations are derived in the body-fixed axis system. But, the position 
of the aircraft must be described by the earth-fixed coordinate system. This is accomplished by 
three consecutive rotations (whose order is important). From Figure 4, the following rotations 
are made: 

1 . Rotate the X\Y\Z\ body coordinate frame about the Z\ axis over the yaw angle ($) to the 
X 2 Y 2 Z 2 body coordinate frame. 

2. Rotate the X 2 Y 2 Z 2 body coordinate frame about the Y 2 axis over the pitch angle (0) to the 
X 3 Y 3 Z 3 body coordinate frame. 

3. Rotate the X^Y^Zj body coordinate frame about the X 3 axis over the roll angle ($) to the 
XY Z inertial coordinate frame. 
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These angles (yaw, pitch and roll) are known as the Euler angles. 

The velocity components between the body-fixed coordinate system and the earth-fixed 
coordinate system are related by a set of orthogonal transformations (a more complete derivation 
may be found in [1] and [2]). These transformations are: 


y 


CqC<% 

S$SqC\$ — C$Sq? 

CJ<£ Sq C\jr S(f> S\$ 

m 


CqS^ 

S$SqS\% + C$C\fr 

— S$C\% 

v 

(18) 

-Sq 

S$Ce 

c$c& j 

[w] 



where C<$ = cos(<$), 5$ = sin('P), etc. 

The kinematic equations relate the Euler angles ($, 0 and $) and the angular velocities (P, 
Q and R). From Figure 4, Co must equal the vector sum of the time rate of change of the Euler 
angles about the k\, j 2 , ami it axes: 


Co = Vk 2 + 0J3 + $i* (19) 

Using transformations similar to Equation (18) and from Equation (13), the angular velocity 
may be written as: 


Co = Pi + Qj + Rk = i 4 f sin 0 + + 

j cos 0 sin $ + 0 cos (20) 

k cos 0 cos $ — 0 sin 

Equating components, the kinematic equations become (in matrix form): 


(p 1 


f 1 

0 

— sin 0 > 


' 4 ' 


Q 

► = < 

0 

cos $ 

cos 0 sin $ 

> < 

0 

> (21) 

UJ 


,0 

— sin <£ 

cos 0 cos <£ 





The time histories of the Euler angles are determined by inverting the 3x3 matrix in Equation 
( 21 ): 


$ = P + Q sin $ tan 0 4- R cos <5 tan 0 

0 = Q cos $ — R sin $ (22) 

4' — (Q sin + R cos $) sec 0 


7 




Figure 4 Rotation from Earth-Fixed to Body-Fixed Frame (from Roskam [1]) 


Aerodynamic Nomenclature 

All forces and moments developed are based on the stability axes system defined by Figure 
5. This system is utilized since experimental data from wind tunnel tests are presented in the 
stability axes. The stability axes are obtained by rotating the body axes ( XY Z ) about the Y = Y, 
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axis. This is done over a rotational angle a (known as the angle of attack) until the body-fixed 
axis (X-axis) coincides with the free stream velocity vector (V pi ). 



The equation for the angle of attack is: 


a = tan 


-l 


W 

U 


(23) 


The equation for the airspeed is: 

V p = ^((7 2 + K 2 + W 2 ) (24) 

The flight path angle ( 7 ) is the difference between the pitch angle and the angle of attack 
(7 = 0 - a). Figure 6 depicts the relationship between the different coordinate axes in the 
longitudinal plane. 
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Figure 6 Relationship Between the Earth, Body and Stability Axis Frame 


Forces and Moments 

To simplify the analysis, all forces and moments are developed in the stability axis system. 
The force ( F 3 ) and moment ( M 3 ) components are: 

F 3 ~ F&* + F 3 aJ* + F% t k a = -Di s + Fyjs - Lk 3 
M 3 = L^i, + Maj s + N^k 5 


where D = Drag, L = Lift and Fy = Sideforce. La, Ma and Na are aerodynamic moments. 

The steady state forces and moments for a straight line flight are assumed to depend only on 
angle of attack (a), sideslip angle (/3), thrust and the control surface deflections of the elevator 
{8e), ailerons (£.4), rudder ( 8 r), and aerodynamic coefficients. These dimensionless coefficients 
are comprised of derivatives evaluated at constant Mach and Reynolds number (e.g.): 


Cd = Cd 0 + Co a oc + Cdi b 8e ( 26 ) 

where: 

• Cd 0 = total airplane drag coefficient for a = 8 e - 0. 

• Coa = total airplane drag change with angle of attack for 8 e = 0. 

• Cq Se = total airplane drag change with elevator angle for a = 0 . 
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A similar analysis is made for the remaining aerodynamic coefficients. 

The aerodynamic forces and moments are expressed in terms of the dimensionless coef- 
ficients, flight dynamic pressure, characteristic length (for the moments only) and a reference 
area: 

D = CoqS = ( 'Cd 0 + Cd„c * + CD, B 8E^qS 

Fy = Cf r qS = ^Cygfl + Cy tA 8 A + Cy fR 8jij qS ^ 

L = CLqS = ( Cl„ + Cz a a + CL tB 8E)qS 


La = CiqSb 

= ^ Ci 0 0 + C lt jA + Cis n s R + Cl P2U~ 3 + ClR 2U~s)^ Sb 
M a = C M qSc 

= ^Cm 0 + C M a a + Cm( b 8e + C Mq — ^jqSc 

Na = CjyqSb 

= (cnJS + Cn, a 8 A + Cn, r 8r + C N p W + CNR-^jj-^jqSb 


( 28 ) 


The flight dynamic pressure is : 



where p is the air density and U, is the aircraft’s airspeed. 
The aircraft thrust vector is also divided into components: 

= T cos (a) 

Ft, =0 

Fj< m = — T sin (a) 

L 3 t = 0 
= —Tdj' 

JVf = 0 


( 29 ) 


( 30 ) 


where d. ? is the moment arm of thrustline. 

For a trimmed airplane, the forces and moments acting on it are in equilibrium. This 
is accomplished when the pitching moment equals zero and when the lift equals the airplane 
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weight. Using Equations (27) and (28): 


0 = (Cm„ + Cm,, « + Cm, b Se^ 

= C Ltrtm = (c Lo + + C LtB 8 E ) 


(31) 


Solving these equations for a and S E determines the trim value for the angle of attack and 
the trim elevator setting: 


Gtrim — 


= 


C L iri , n C M s b ~ c m, b + C L(b C M ft 

CL a C M(B - C L(b C m „ 

C L ^C Mn -CL 0 C Ma + Cl„ Cm 0 


Cl>C m „ - C La C M , K 


(32) 


(33) 


Combining the rigid body equations and the aircraft kinematic equations yields a full set 
of flight dynamics equations that describe any aircraft flight path or motion. To accomplish 
this, all forces must first be transformed from the stability axis system to the body axis system. 
From Figure 5: 


F b = —D cos (a) + L sin (a) + T cos (a) 

F b y = F '• (34) 

F b z = — D sin (a) — L cos (a) + T sin (a) 

The force equations (14) are then solved in terms of the linear velocity rates (U, V and 
W). Also, the moment equations (17) are solved in terms of the angular velocity rates ( P , Q 
and R ). To simplify the analysis, let: 

k i = Ixz{Iyy ~ Ixx ) 

^2 = hz(hz ~ Wy) 
kz = Ixzihz - Wy) 
k* = IxxVyy - Ixx) 

The full flight dynamic motions are now expressed as a 12 th order set of nonlinear differential 
equations (three of which are simple integrations of the linear velocities to yield linear positions). 
These equations can now be integrated to yield flight trajectories, after a transformation is made 
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to the earth-fixed coordinate system by Equation (18), of any variable (e.g. 
airspeed, etc). 

plane altitude, 

U = —g sin (0) H — + V R — WQ 

m 

(36) 

F b 

V = g sin ($) cos (0) -1 — — UR + W P 

m 

(37) 

F b 

W = <7Cos($) cos(0) -1 — — + UQ — VP 

m 

(38) 

^ IxzNa + IzzLa - hPQ - l\zQ R + txzhzPQ - fc 2 RQ 

Ixxhz - Ijcz 

(39) 

. M a - {Ixx ~ hz)PQ ~ Ixz(P 2 ~ r2 ) 

Iyy 

(40) 

• IxxN a + IxzLa + I 2 xz p Q ~ RQ - *4 PQ - IxxIxzQR 

Izzlxx - I\z 

(41) 

Q — Q cos $ — R sin $ 

(42) 

= (Q sin $ + R cos $) sec 0 

(43) 

$ = P + Q sin $ tan 0 + R cos $ tan 0 

(44) 


Table 1 Flight Dynamic Equations 


The full flight dynamic equations in Table 1 may be linearized into second order differential 
equations by utilizing small disturbance theory. This theory applies to small deviations (for 
angle of attack, sideslip and control surface deflections, etc.) relative to some steady state flight 
condition. For an automatic landing system, this assumption is valid since only small angle 
deviations caused by turbulence or control variables are encountered. This theory is also useful 
in analyzing the stability of an autopilot by using longitudinal transfer functions. 
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Small Perturbation Equations (stability axis) 

Perturbed state equations are derived by replacing all motion variables by a steady state 
value and a perturbation: 


U = U 0 + AU 
P = P 0 + AP 
$ = #0 + A$ 


V =V 0 + AV 
Q — Qo + AQ 
0 = 0 O + A0 


W = W 0 + AW 
R — Rq -F A R 
$ = + A<t> 


This also applies to all forces and moments: 

F x = F Xo + A F t F y = F yo + AF y F z = F Zn + A F z 

M = M 0 + AM L = Lq + AL N = N q + AN 


(45) 


(46) 


The change in forces and moments can be expressed in terms of the perturbation variables 
by using a Taylor series expansion: 


dF : 


8F : 


dF x 


8F X 


A F x = — f A U + ^AW + — ± A 6 S + ^A5 T 


dU 


dW 


88 E 


86? 
dF 




AF Z = ^A U + ^AW + ^A Q + ^-A S E + ^A 8 T 


dU ' dW 


. , . w Arr dM 
AM “ dU A + 8W 


dQ 


AW + ^AQ + 


dQ 


dS E 88? 


8M 8M 

—~A8 e + -r-r-A 8 t 
(JOE OvT 


(47) 


8L 8L 8L 8L A r 8L A r 

A L = j-AV + —A P + — A R + —A 8 R + — A£ 4 
8V 8P 8R 88 e 88 a 


_ 8N A T , 8N A „ 8N A „ 8N A _ 3iV A . 

AN= _ AV+ + —A* + —ASn + 

The partial derivatives in Equation (47) are called the stability derivatives. For convenience, 
these derivatives are divided by the aircraft mass. These new symbols are defined as dimensional 
stability derivatives (e.g.): 




( 48 ) 
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Replacing the forces and moments with Equation (25) the longitudinal dimensional stability 


derivatives, in the stability axis 

system, become: 



-2C Dn qoS 

mUo 

-2C D(m q 0 S 

Xs ■- 

-2CL a qoS 
u ~ mU 0 


r? -CL h qQSc 

“ 2 mUo 

„ -Ci Q qoSc 

q ~ 2 mUo 

„ -Cl,JoS 

Z ‘’ = m 


2CM„qoSc 
u ~ IyyUq 

,, -CM„qoSc 

M a = 

lyy 

,, CM„qoSc 2 

Mi= 2 IyyU 0 

(49) 

Cmo^oSc 2 

q ~ 2 I yy Uo 

,, qoSc 

M 6b = t b 

lyy 

,, Cm t JoSc 

M T a = j 

lyy 


Y ~{Cd* - 

vV a — 

- C Lo )qoS v 

~{Ci a + Con'jqoS 



These dimensional stability derivatives, along with the force and moment perturbation 
equations, can now be evaluated into the full flight dynamic equations (Table 1). To simplify 
the analysis, for small angle deflections, let: 


cos 0 ~ 1.0 
sin © ss © 


(50) 


Taking a Laplace Transformation of the resulting equations determines the longitudinal 
transfer functions with the elevator setting (6e) as the input and the horizontal velocity component 
( U ), angle of attack (a) and pitch angle (©) as the output variables. These longitudinal transfer 
functions, in matrix form, are: 


Mi M2 Mz ) ( 6^(3) 
Mi ^22 ^23 > l > 

-431 -432 -433 J l ^ > 


Xf. ) 
Z 6b 

m 6b J 


(51) 
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where: 


An = (s- Xu) 

An = -X* 

vli 3 = g cos @o 
A 21 = —Z u 

A 22 = s{U Q - Z 6 ) - Z a (52) 

A 23 — — (Zq + Uo)s -f- g sin 0o 
An = —(M u + Mt„) 

^32 = ~(sMa + M a ) 

An = s 2 - sM q 


Taking the inverse of Equation (51) yields the transfer functions for the longitudinal mode. 
Using Cramer’s Rule, the denominator of the transfer functions is a fourth order polynomial. 
The roots of this polynomial form the characteristic equation and determine the stability of an 
aircraft. The characteristic equation is represented by two oscillatory modes of motion (short 
period and phugoid). 


Short and Long Period Approximations 

The longitudinal motion of an aircraft is determined by two modes. The first is the short 
period mode. This is characterized by a highly damped, high natural frequency oscillation at 
approximately constant speed (AU « 0). The other mode is the long period or phugoid mode. 
This is caused by a gradual change between potential and kinetic energy and is characterized 
by a low damped, low natural frequency oscillation at approximately constant angle of attack 
(Aa a; 0). Both modes are determined by the characteristic roots of Equation (51). The concept 
of these modes is illustrated in Figure 7. 
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The short period mode is obtained by approximately neglecting the velocity term in Equation 
(51). This equation now reduces to: 


sUo-Z a -sUo 

-Ma-M a s 2 -sMJX^S \M 6b J 


(53) 


Applying the inverse to this equation: 

ajf) _ sZ 6b + (Ms b Uq - MgZss) 

6e{s) ~ s{s 2 Uq - s(M q Uo + Z a + UoMa) + ( Z a M q - M a U 0 )} 

Q(j) 3 (Uq^Sb + Zs^Mg) + (MoZsb ~ Z a Ms E ) 

S E (s) ~ s{s*U 0 - s{M q U 0 + Z a + U q M 6 ) + ( Z a M q - M a U 0 )} 


(54) 


Comparing the denominator of this equation to the standard frequency form of ( s 2 + 
s2(ujnsr + ^nsp)’ the natural frequency and damping ratio of the short period mode are: 
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u ’nsr 



- M a 


Csp 


z c M£o + Z a + MM) 

2 Wnsp t/0 


( 55 ) 


The phugoid approximation is derived by neglecting the angle of attack term in Equation 
(51). The resulting equations for this mode become: 


Ujs) _ Zg.Upg 

6 e{ s) (s 2 U 0 — sX u Uq — Z u g ) 

0 ( 3 ) = -Z Sb U*{s-X u ) 

Se{s) ( s 2 Uo - sX u Uo - Z u g) 


(56) 


The natural frequency and damping ratio for the phugoid mode are: 





Cp 


-X u 

2 U np 


(57) 


The combination of the short period and phugoid modes describes the complete longitudinal 
aircraft motion (i.e.): 


0 ( 3 ) ^©^(70! s + i)(r e ,3 + 1 ) 



(58) 
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